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Abstract. We consider the paths of a Gaussian random process x(t), x(0) = not 
exceeding a fixed positive level over a large time interval (0,T), T> 1. The probability 
p(T) of such event is frequently a regularly varying function at oo with exponent 0. In 
applications this parameter can provide information on fractal properties of processes that 
are subordinate to x(-). For this reason the estimation of 8 is an important theoretical 
problem. Here, we consider the process x(t) whose derivative is fractional Brownian motion 
with self-similarity parameter < H < 1. For this case we produce new computational 
evidence in favor of the relations logp(T) = — 6*logT(l + o(l)) and 9 = H(l — H). The 
estimates of are to within 0.01 in the range 0.1 < H < 0.9. An analytical result for the 
problem in hand is known for the markovian case alone, i.e., for H — 1/2. We point out 
other statistics of x(t) whose small values have probabilities of the same order as p(T) in 
the log scale. 
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1 Introduction 



The asymptotics of tail probabilities P(maxx(s) > u), u — > oo of a Gaussian 

random process x(t) is a problem that has received a sufficiently complete 
solution (see, e.g., [12] )• Considerable progress also occurs in the study of 
probabilities of small deviations, that is, of events of the form {max \x(s)\ < 

e} as e j (see, e.g., [7j, [H]). However, there is practically a complete 
absence of general approaches to the analysis of probabilities of unilateral 
small deviations, i.e., of probabilities of the form 

p{u\T) = P(max(x(s) - x(0)) < u) (1) 

for small u or large T. Here, T is the similarity ratio and A a finite closed 
interval that contains 0. If x(-) is a self-similar process, that is, finite dimen- 
sional distributions of {x(tA)} and {lA^a^t)} are identical for any A / 0, 
then p(u\T) = p(uT~ h \l), so that the problems on the asymptotics of (0) for 
the u and T indicated above are equivalent. 

Exponential asymptotics are typical of large and small deviations of x(-), 
while power law ones are typical of unilateral small deviations. Power law 
asymptotics are rather popular in physics, since they frequently provide in- 
formation on fractal properties of physical processes (see the examples later 
in this paper). In this connection one is faced with the task of calculating 
the exponent 

6 = -limlogp(l|T)/logT, T-+oo. (2) 

We shall refer to the exponent 9 for a sequence of events At in what 
follows, when log P(At) = — 01ogT(l + o(l)) as T — ► oo. 

There are few explicit estimates of (J2J). Molchan (PJ) has found the ex- 
ponent 9 for fractional Brownian motion (FBM), bn{s). That one-parameter 
family of processes is specified by the requirement bn(0) = and by the 
structure function 

E\b H (t)-b H (s)\ 2 = a\t-s\ 2H , 0<H<1, 
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where H is the self-similarity parameter of the process. When H = 1/2, 
bn{s) becomes Brownian motion. 

According to 0, 9 = 1 — H when A • T = (0, T) and 9 = 1 when 
A-T = (— T, T). The last estimate is of interest in that it is independent of H. 
The property in question also remains valid for FBM with multidimensional 
time, t G R d . In that case one has 9 = d, if A • T = {t : \t\ < T}. 

In mathematical physics one is interested in the parameter 9 for the in- 
tegral of fractional Brownian motion (IFBM), that is, for the process x(t) = 
fo^H{s)ds, \t\ < oo [H], [TH]. By now, the problem of 9 has been solved for 

H = 1/2 only, when {x(t), x (t)) is a Markov process. Sinai ^3] showed 
that 9 = 1/4 for A ■ T = (0, T). Since paths of x(t) in (-T, 0) and (0, T) are 
independent, it follows that = 1/2 for A-T = (— T, T). The general case 
< H < 1 was studied by the present authors ([101) both analytically and 
numerically. Our analysis suggests the following hypothesis: 



9=r~ H > A-^ = (-T,T) (a) 

A-T = (0,T) (6). U 

The hypothesis (J2t) is corroborated by a related result in [3] and by numerical 
calculations. As to (J2b), it has been confirmed numerically only for the 
interval 0.1 < H < 0.6 to a low accuracy (5 = 0.03). Here we continue 
the numerical analysis of our hypothesis (JHt 1 ) by making use of analytical 
results derived in [TO], but radically modify the evaluation strategy. This 
allows corroboration of (3b) in the entire range of H to within 5 < 0.01. One 
by-product is to provide support in favor of the following asymptotics: 



P(x(s)x(l) > 0, 1 < s < T) = 0((\ogT)^ H) T- H ^- H) ), T -> oo, 
where a(H) may have the form H — 1/2. 



2 Examples 



Consider a few examples where exponents like (2) are used. 

Let M(t) = max(x(s) — x(0)) be the record function of the process x(-). 

Levy has shown (see 0]) that the record function of Brownian motion is 
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similar to Cantor's staircase. Its points of growth make a set S of Hausdorff 
dimension dimS 1 = 1/2. The situation for the FBM family is analogous: 
dimS* = H. This result is based on two fact for FBM in A • T = (0, T): the 
exponent is 9 — 1 — H , and 



log[P(G(T) < 1)/P(M(T) < 1)] = o(logT), (4) 

where G(T) is the position of the maximum M(T) in the interval A • T ( 9J). 

Denote by Z + the first zero of bu{t) after the time t = 1; then, similarly 
to (H, 



\og[P(Z + > T)/P(M{T) < 1)] = o(logT). 

Consequently, the log asymptotics of P(M(T) < 1) also determines the log 
asymptotics of long excursions of FBM. The interest in such asymptotics is 
rather broad in physics and engineering (see review [S]). 

There is a similar problem in geology bearing on the dimension of fractal 
sets. Molchan and Turcotte ([11) consider a simple one-dimensional model 
of sedimentation in shallow seas. The model involves two mechanisms that 
produce the sediment: tectonic forces which cause relative rise/fall of sea level 
on the one hand and erosion on the other. A sea level fall brings the upper 
sedimentary layers above the water and causes a fast (on the geological time 
scale) erosion of these. The erosion process causes time gaps (unconformities) 
in datings of the sedimentary layers. A sea level rise provokes sedimentation. 
The latter process is considered to be rapid, so that the sediment follows sea 
level variations practically continuously. Supposing the sea level history to 
be described by a process x(t), the drilling data (layer depth and date of 
formation) must be described by the function 



y(t) = minx(s). 

s>t 

Points of growth of y(t) (the set S) correspond to the dates of those layers 
which have been preserved in the sedimentary sequence during the entire 
history of sedimentation. The common assumption makes x(s) a process 
with stationary increments; the self-similarity of S in an extensive range of 
scales suggests self-similarity for the stochastic component of x(-). When x(s) 
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is also assumed to be Gaussian, one arrives at the model x(s) = as + bn{s), 
a > (see e.g. j2]). The dimension of the support of y(t) for the above 
model of x(-) can be found similarly to the preceding problem, giving the 
result dim 5 = H (PU). 

The second example is concerned with the fractal nature of solutions to 
the Burgers equation with random initial data: 

u t + uu x = vu xx (5) 
u(Q,x) = v(x), \x\ < oo 

where J v(a) da = o(x 2 ), i>1 (see, e.g., [TE]). 

When the viscosity v is infinitely small, equation (J3J) when considered 
in the one-dimensional case describes the following dynamics of adhesive 
particles. A particle at the initial time at position a has the mass da and the 
momentum v(a)da. The particle starts moving at velocity v(a) and preserves 
that velocity until the first collision with its neighbors. The colliding particles 
stick together and continue their motion by following the laws of preservation 
of mass and momentum, and so on. The question that arises is what is the 
dimension of the initial positions of the particles S that have not collided 
until time to- Such positions are called Lagrangian regular points. 

The question is amenable to a purely geometrical interpretation. Let U(x) 
be a convex minorant of the curve £(x) = J? v (a) da + x 2 /2. Then the right- 
hand derivative U'(x) is nondecreasing function, and its points of growth 
correspond to Lagrangian regular points S. The dimensionality problem of 
S in the case v(a) = 61/2(0) has been solved by Sinai [TJ]: dim S = 1/2. The 
solution is based on estimating the exponent (2) for the integral of Brownian 
motion in the interval (0, T). For the case v(a) = b H (a) it can be asserted 
that dim S = H, if B for IFBM in A ■ T = (— T, T) admits of the bound 

0>{\- II) mr 

The exponents (J2J) for the maximum of IFBM are closely related to the 
exponents of other statistics. To be more precise we introduce the following 
notation: M(T) = maxx(t), G(T) is the position of M(T) in A • T, and 

A(T) = J la;( s )>o ds is the occupation time of x(t) above in A • T. 

A-T 

Statement 1. (Molchan & Khokhlov, [JU]). Let A ■ T = (0,T) or 
(-T, T). For x(t) = IFBM(t), t G A-T the exponents of the events {MiT) < 
1}, {\G(T)\ < 1}, {AiT) < 1, \G{T)\ < T}, {x(t) < 0, t G A • T, |t| > 1}, 
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(a) are identical, when they exist; 

(b) exist or do not exist simultaneously. 

It is a known fact that G(l), A(l) and the last zero, Z, for Brownian 
motion in (0,1) have identical distributions, namely, the arcsine law. State- 
ment 1 is also true for FBM ([H]). Therefore we can consider this result as a 
weak version of the arcsine law for FBM and IFBM. Statement 1 provides a 
certain degree of freedom in numerical analyses of (J2J). 

3 Evaluation of 6 

We are going to evaluate the exponent of the event {M (T) < 1} for the IFBM 
process in the interval (0,T),T>1 using the Monte Carlo method. To select 
a suitable strategy note the following. Statement 1 gives some information 
on the possible structure of a typical IFBM path with a low maximum, 
M(T) < 1. Namely, the position of the maximum G{T) and the total time 
A(T) where IFBM > do not practically grow (the growth is most likely to 
be a logarithmic one). In most of the cases the path goes to the lower half- 
plane after the lapse of a fixed time, because there are no limitations on the 
amplitude from below. Consequently, the essential information concerning 
the low maximum is available around the initial point t = 0. 

Suppose we generate IFBM on a uniform lattice (0,5, ... , L5) using tri- 
angular factorization of the correlation matrix. We note that the correlation 
structure of the process should be reproduced exactly in the case under con- 
sideration, since we are dealing with rare events {M(T) < 1}. Hence we 
shall need a memory of order L 2 for generating a large number of samples. 
However, the second-degree growth in memory leads to limitations on lattice 
width, hence on the information concerning the low maximum. This heuris- 
tic argument is borne out in JU] . Following the above strategy, we have not 
succeeded in evaluating 9 for 0.6 < H < 1. 

The way out consists in considering the IFBM in a log time scale. To do 
this, we consider the Lamperti transformation which converts {IFBM(t),t > 
0} into a stationary process x(t): 

x(t) = cexp(-(l + .H»rFBM(e*), \t\ < oo. (6) 

We have normalized the process so as to make Ex 2 (t) = 1, hence c = 
(1 + 2H)~ 1 / 2 . The correlation function r(t) for x(t) is 
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r(t) = 2c 2 (l + H) cosh(Ht) - c 2 cosh((l + H)t) + 0.5c 2 1 sinh(t/2)| 2 " +2 
and 



r{t) 



1 - 0.5(1 - H 2 )t 2 + 0.5c 2 t 2+2H + 0(t 4 ), t -> 
c(H) -exp(-pi) • (l + o(l)), t^oo, 



where p = min(if, 1 — iJ). 

We can see that has a smoothness of order 1 + H — e (e > 0), the 
number of zeroes of x(-) is locally finite, and the mean interzero distance is, 
according to Rice, given by 

A = n(l-H 2 )- 1 / 2 . (7) 

The transformation © does not preserve the point (G(T), M(T)) as an 
extreme one in a sample of x(-). However, the event 

{IFBM(t) < 0, 1< t < T} = A T 
is easily rewritten in terms of x(-) on a finite interval: 

^l T = {x(r)<0, < t < T' = InT} := A T >. 
Consequently, the desired exponent 



6 = - lim hi P(A T ')/T' (8) 

is converted to exponential from a power- law one. As a result, the evaluation 
of 9 splits into three steps: 

• generating a stationary process x(t), t > 0; 

• evaluating the distribution of the first zero, Z, for x(t); 

• finding (jHJ) for the tail of the distribution of Z. 
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The Generation of x(t). The process x(t) was generated as a stationary 
sequence x(5 ■ k) with the exact correlation function r(5k), k = 0, 1,...,L. 
This was done using the triangular representation 



x(kS) = a(i\k)si 



(9) 



t=0 



in terms of the standard white noise The representation is implemented 
by using the progressive Schur algorithm (see [I]). The discretization step S 
is specified by the number of points n per mean period A for the zeroes 
of x{t) (see flJJ)). Since x(-) has a smoothness of order ~ 1 + H, and we are 
interested in long excursions of x(t), we can well use moderate numbers for 
no; we had no = 50 in the calculations. The length L can be found from the 
requirement P(Z > L5) = e where e is small. 
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Figure 1: Empirical functions y — — In P(Z > x) based on N = 100,000 
paths of IFBM: H = 0.1 0.4 (left) and H = 0.5 ^ 0.9 (right). Straight lines 
correspond to y = H(l — H)x. 

As is shown by some preliminary evaluations of the distribution of Z (see 
Fig. HJ, the function P{Z > t) for H = 0.1 0.9 is well fitted with an ex- 
ponential function: exp(— H(l — H)i). Consequently, one has approximately 
L5 ~ Z e , where 



Z £ = - In £/[#(!-#)]. 



(10) 
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When e = 1(T 4 and H = 0.2-0.8, one has L5 ~ 60, and L5 = 100 when H = 
0.1, 0.9. The Schur algorithm is a recursive one, so it may become unstable 
as L increases. The instability manifests itself in parasitic oscillations of the 
a(i\k) at large i (see ©). These effects are typical of i-5 > 40 and hardly can 
always be overcome by using available accuracy (for instance long double 
in C) . Hence it follows that 

• there are computational difficulties in the way of analyzing the distribu- 
tion of Z when H is close to either or 1; more exactly, when \H — 0.5 1 > 0.4; 

• evaluation of 9 from the values Z > Z £ , e < 10~ 4 calls for higher 
computation accuracy. 

The last conclusion is important, because we do not know when the log 
linear asymptotics for the tail of the distribution of Z becomes valid. 

Evaluation of 9. The parameter 9 was evaluated in the series of intervals 
A e = (Z e ,Z e /va) with e = 0.01; 0.003; 0.001 for the range H — 0.1 — 0.9. 
Figure ^ does not contradict the assumption of linearity for the plot of 
(t, In P(Z > t)), t G A s . For this reason we use for the slope of the plot 
the maximum likelihood (ML) estimate 9 corresponding to the distribution 
{c(9)e~ 0t ,t G A e }. Namely, 9 = x/\A £ \ where |A e | = j- In 10 is the length of 
A e , #o = H(l — H), and x is the root of 

x- 1 - (e* - l)- 1 = [<Z> £ -Z £ ]/\A £ \. (11) 

Here, < Z > £ is the empirical mean of all Z observed in the interval A e . 

ML estimates of 9 based on N = 16 x 300, 000 paths of x(-) are shown 
in Fig. 121 They demonstrate that the slopes 9 £ in the intervals A £ , e = 
0.01,0.003,0.001 are well consistent among themselves and are identical with 
the hypothetical values 9q = H(l — H) to within 0.01. Note that the left- 
hand endpoint of A e is approximately identical with the (1 — e) quantile of 
the distribution of Z. Consequently, the number of observations N £ used to 
estimate 9 £ is approximately equal to N ■ e. 

The ML estimate of 9 for the truncated exponential distribution 
{c{9) exp(— 9t), t G (hie -1 , lnlO/e)^ -1 } has the standard deviation a £ ~ 

— 1/2 

1.7014 9N £ , where N £ ^> 1 is the number of observations. In our case the 
slope estimates are close to #o — H(l — H), while N £ is large, hence a £ = 

— 1/2 

1.7 9qN £ can serve as a satisfactory theoretical estimate of the standard 
deviation for 9 £ . This statement is corroborated by our experiments carried 
out to check the operation of the random digit generator and the Schur 
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Figure 2: Hypothetical 6q = H(l — H) and empirical 9 exponents for three 
intervals A e , e = 0.01, 0.003, 0.001. Number of paths of IFBM is iV = 
16 x 300,000. 



algorithm. The estimates 9 £ were derived above by averaging over 16 serial 
estimates 9f v . Each series consists of 300,000 paths of x(-). The empirical 
variance of the averaged estimate 9 £ is in good agreement with the theoretical 
value. Examples are given below in Tabled 

Figure El shows interval estimates 9 £ ±a £ of the slope 9 e for e = 0.01. Even 
though we have seen above that \9 £ — 9 \ < 0.01, Fig. El provides evidence 
of a significant discrepancy between empirical and hipothetical estimates 
of 9. Furthermore we can see that 9 £ (H) > 9 (H) when H < 0.5 and 
9 £ {H) < 6 (H) when H > 0.5 for all e = 0.01, 0.003 and 0.001 (Fig. EJ- 

The above inference cannot be ascribed to the effect of discretization. 
This is confirmed by the estimates of 9 for H = 0.7 with different values of the 
discretization parameter: no = 30 and 100 steps per period Ao (see Table [Q. 
The estimates of 9 £ correspond to the intervals A e = (In-, In—) ■ 9^ with 
£ = 0.1, 0.01, 0.001 and were derived by averaging over r serial estimates 
#^ er . Each series consists of 10 6 paths of x(-), r = 120 in the case n = 30 
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n 


e 


% 


a 9 




.100 


.2002 


.00009 


.00011 


30 


.010 


.1998 


.00031 


.00032 




.001 


.1997 


.00089 


.00096 




.100 


.2001 


.00011 


.00012 


100 


.010 


.2001 


.00040 


.00037 




.001 


.1989 


.00107 


.00112 



Table 1: Estimates 9 in intervals A e for H = 0.7 and discretization parameter 
n = 30 (top) and no = 100 (bottom). The standard deviations of 9 are erg- 
(empirical) and <r~ (theoretical). 

and r = 90 in the case n = 100. Tabled also lists empirical and theoretical 
standard deviations of the estimates of 9 £ . It appears from Tabled that the 
estimates 9 £ « 0.200 are independent of the discretization step and there are 
significant deviation of 9 £ from the hypothetical value 9q = 0.21. 

We thus have to reject the hypothesis that the tail of the distribution of 
Z has a purely exponential asymptotics with parameter 9 = 9 . However, 
this does not rule out the hypothesis proper of the exponent 9 = 9 for the 
events {Z > T}, T 3> 1. We are going to show that the distribution 

P(Z > t) = ct a e- 9ot (l + o(l)), t^oo (12) 

is consistent with our estimates of 9 £ . 

To do this, let us replace the empirical mean < Z > £ in (jllj) with EZl ze /^ s 
assuming (fT2j) and o(l) = 0. Solving the equation yields the expected value 
of 9 £ under (JT2J). These estimates 9 £ are shown in Fig. El for a = H — 0.5. It is 
seen that the 9 £ are in very good agreement with the empirical estimates 9 £ 
for all H = 0.1 — 0.9. Figure 0] provides a more detailed view of the residuals 
R £ = 9 £ — 9 £ . ^ 

Both Figs. 3 and 4 show that the empirical estimates 9 £ can be well fitted 
using the extra parameter a. However, it is very difficult to get a with a 
suitable resulution 5, say 5 = 0.02 — 0.03, if a is small, as is the case for the 
model a = H — 0.5. Indeed, suppose the parameter 9 is known and 9 = 9 Q . 
The Cramer-Rao inequality yields the optimal variance t for a: 
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H 

Figure 3: Hypothetical values of 6q = H(l — H) and interval estimates 
9 £ ± a £ of the exponent 9 for the interval A e , e = 0.01. Number of paths is 
N = 25 x 300, 000. Boxes correspond to the expected values 6 £ of 8 under 
condition H2*l and a = H — 0.5. 



a opt ~N^ /2 a- l (\nZ A ) 

where a 2 (lnZ A ) is the variance of InZ for observations of Z in the interval 
A, Nj\ being the number of the observations. Hence 

N A ~[a opt -a(\nZ A )]- 2 . 

If A = (Z e , A/6 ) and Z £ is given by (fTU|). then the total number of paths of 
x(-) is 

N^[a opt -a(\nZ A )]- 2 s-\ 
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0.005 r- 
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-0.003 - 

-0.004 - 

-0.005 I i i i ' i i i i i 

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 

H 

Figure 4: Residuals R = 9 e — 9 £ , e = 0.01 (see Fig. EJ). Vertical lines 
correspond to two levels of 9 £ deviations: a)±cr e (bold) and b)±2d £ , (thin) 
where a £ is empirical standard deviation of the estimate 9 

One has |a| < 0.5 in the model a = H — 0.5. For this range of a with e = 0.01 
and 0.001 for A = 20 one has a(lnZA) — .15883 and 0.11536, respectively. 
The requirement <T op t = 0.01 makes the number N large enough: 40 • 10 6 if 
e = 0.01 and 750 ■ 10 6 if e = 0.001. 



4 Conclusion 

We have shown that the exponent hypothesis #o = if (1 — H) for a series 
of statistics related to IFBM (see Statement 1) is well corroborated by our 
computations. Further refinement of the tail probabilities for these statistics 
faces considerable computational difficulties in view of the amount of compu- 
tation required, the computation accuracy, and checks on the random digit 
generator. 

The symmetry of the exponent: 9(H) = 9(1 — H) which 9 has by defi- 
nition is not obvious in our problem because of quite different properties of 
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FBM for H < 1/2 and H > 1/2. Our analysis suggests that this difference 
can manifest itself in more sophisticated asymptotics, e.g. 

P(IFBM(t) < 1, < t < T) = 0(T- e °(\og T) a{H) ) 

where a(l/2) = 0. Of course this relation, as well as 6q = H(l — H), are 
needed in analytical corroboration. 
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